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Abstract 

We are concerned with the numerical solution of a class integro- 
differential equations, known as Neural Field Equations, which describe 
the large-scale dynamics of spatially structured networks of neurons. 
These equations have many applications in Neuroscience and Robotics. 
We describe a numerical method for the approximation of solutions 
in the two-dimensional case, including a time-dependent delay in the 
integrand function. Compared with known algorithms for this type 
of equation we propose a scheme with higher accuracy in the time 
discretisation. Since computational efficiency is a key issue in this type 
of calculations, we use a new method for reducing the complexity of the 
algorithm. The convergence issues are discussed in detail and a number 
of numerical examples is presented, which illustrate the performance 
of the method. 


1 Introduction 

In recent years significant progress has been made in understanding the brain 
electrodynamics using mathematical techniques. Neural field models repre¬ 
sent the large-scale dynamics of spatially structured networks of neurons in 
terms of nonlinear integro-differential equations. Such models are becoming 
increasingly important for the interpretation of experimental data, including 
those obtained from EEG, fMRI and optical imaging [5]. These equations 
also play an important role in Cognitive Robotics, since the architecture of 
autonomous robots, able to interact with other agents in solving a mutual 
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task, is strongly inspired by the processing principles and the neuronal cir¬ 
cuitry in the primate brain (see [7]). All this explains why Neural Field 
Equations in several dimensions are a very actual and important subject of 
research. 

Moreover, simulations play a fundamental role in studying brain dy¬ 
namics in Computational Neuroscience, and to understand diseases such as 
Parkinson, as well as the effect of treatments, such as in Deep Brain Stim¬ 
ulations (DBS) or Transcranial Magnetic Stimulations (TMS). Thus, the 
availability of efficient, fast, reliable numerical methods is an important in¬ 
gredient for improving the effectiveness of techniques such as DBS or TMS in 
many therapeutic applications. Integro-differential equations in several spa- 
cial dimensions are quite a challenge for numerical simulation, because the 
standard approaches require a very high computational effort. Maybe this 
explains why there exist very few studies concerning the numerical analysis 
of NFEs. Some important papers which inspired our work are Eicaci]; 
as other relevant references we can cite mm . Overall the lack of efficient 
algorithms represents a serious drawback for the use of NFEs in practical 
applications. This is the main motivation for the present work. 

We are concerned with the numerical solution of the following integro- 
differential equation: 

c^D(x,t) = - V{x,t) + j^K{\x - y\)S{V{y,t))dy, (1) 

t e [0,r],x G D C 

where the unknown V{x,t) is a continuous function F : D x [0, T] —?• M, I, 
K and S are given functions; c is a constant. In this article, by \x — y\ we 
mean ||a: — y\\ 2 - We search for a solution V of this equation which satisfies 
the initial condition 


F(x,0) = Fo(^), x£n. (2) 

Along with equation Q we will also consider 

c^I^(®,i) = - V{x,t) + j^K{\x - y\)S{V{y,t- T{x,y))dy, (3) 

t G [0,T], X G D C M^, 

where r(x, y) > 0 is a delay, depending on the spatial variables. In the latter 
case, the initial condition has the form 

V{x,t) = Vo{x,t), X G D, t € [-Tjnax,0], (4) 
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where Tmax = maxs,i 7 en 'r{x, y). By integrating both sides of Q with respect 
to time on [0,T], we obtain the Volterra-Fredholm integral equation: 


cV{x,t) = Vo{x)+f ({I{x,s)-V{x,s)+[ K{\x-y\)S{V{y,s-T{x,y))dy 
Jo V Jn 

(5) 

t e [o,r],x G n c 


The existence and uniqueness of a solution of equation 0 in the case 
Q = m > 2, was proved in m, both in the case of a smooth and dis¬ 
continuous function S. An analytical study of equation Q was carried out 
in [8], where the authors have addressed the problems of existence, unique¬ 
ness and stability of solutions. They define the Banach space C = C{[l,F), 
were I is some real interval, containing 0, and F denotes the Banach space 
L^(n,M),with the norm 


ds, 


|'h||p’ = 



'I'2(r)dr, 


VT G F. 


Hence the norm in C is defined by 


<1>||^ = sup. 


t£l 


<h2(r, t)dr, 


G C. 


They have proved that if iF G I G C{[—Tmax , cx)[, T) and r G 

(7(0^, ]R_i_), then for every Vq £ C{[—Xmax,0]) equation Q has a unique 
solution V G oo[, F) n C{[—Tmax,oo[,F). Subsequently in this paper 

we will use the same notations and definitions of norms. Moreover, as the 
authors of [8], we will assume that S and its derivative S' are positive and 
bounded. When solving numerically equations of the forms 0 and ([^, they 
are often reduced to the form ^ ; therefore we begin by discussing literature 
on computational methods for Volterra-Fredholm equations. Starting with 
the one-dimensional case, without delay, Brunner has analysed the conver¬ 
gence of collocation methods [3], while Kauthen has proposed continuous 
time collocation methods m- In [9] an asymptotic error expansion for the 
Nystrom method was proposed, which enabled the use of extrapolation al¬ 
gorithms to accelerate the convergence of the method. Another approach 
was developped by Z. Jackiewicz and co-authors m,m. who have applied 
Gaussian quadrature rules and interpolation to approximate the solution of 
integro-differential equations modelling neural networks, which are similar 
to equation Q. 
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In all the above mentioned papers the authors were concerned with the 
one-dimensional case. In the two-dimensional case, the required computa¬ 
tional effort to solve equations Q and Q grows very fast as the discretiza¬ 
tion step is reduced, and therefore special attention has to be paid to the 
creation of effective methods. This can be achieved by means of low-rank 
methods, as those discussed in when the kernel is approximated by 
polynomial interpolation, which enables a significant reduction of the di¬ 
mensions of the matrices. In [3|, the authors use an iterative method to 
solve linear systems of equations which takes into account the special form 
of the matrix to introduce parallel computation. Concerning equation Q, 
besides the existence and stability of solution, numerical approximations 
were obtained in |8]. The computational method applies quadrature rule 
in space to reduce the problem to a system of delay differential equations, 
which is then solved by a standard algorithm for this kind of equations. A 
more efficient approach was recently proposed in mm. where the authors 
introduce a new approach to deal with the convolution kernel of the equation 
and use Fast Fourier Transforms to reduce significantly the computational 
effort required by numerical integration. 

The above mentioned equations are known as Neural Field Equations 
(NFE) and have played an important role in mathematical neuroscience for 
a long time. Equation Q was introduced first by Wilson and Cowan [T6] . 
and then by Amari [T], to describe excitatory and inhibitory interactions in 
populations of neurons. While in other mathematical models of neuronal 
interactions the function V (membrane potential) depends only on time, in 
the case of NEE it is a function of time and space. The function I represents 
external sources of excitation and S describes the dependence between the 
firing rate of the neurons and their membrane potential. It can be either 
a smooth function (typically of sigmoidal type) or a Heaviside function. 
The kernel function K(\x — y\) gives the connectivity between neurons in 
positions x and y. By writing the arguments of the function in this form 
we mean that we consider the connectivity homogeneous, that is, it depends 
only on the distance between neurons, and not on their specific location. 

According to many authors, realistic models of neural fields must take 
into account that the propagation speed of neuronal interactions is finite, 
which leads to NEE with delays of the form ([^. 

In the present paper we propose a new numerical approach to the Neu¬ 
ral Eield Equation, in the forms 0 and 0. One remarkable feature of our 
method is that we use use a implicit second order scheme for the time dis¬ 
cretisation, which improves its accuracy and stability, when compared with 
the available algorithms. Moreover, to reduce the computational complex- 
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ity of our method we use an interpolation procedure which allows a drastic 
reduction of matrix dimensions, without a significant loss of accuracy. This 
improves the efficiency of the algorithm. 

The outline of the method is as follows. In Sec. 2, we describe the nu¬ 
merical algorithm, for equation Q; its stability, convergence and complexity 
are analysed. In the same section, we introduce an algorithm for equation 
In Sec. 3 a set of numerical examples is presented, which illustrate 
both cases, and the numerical results are discussed. We finish with some 
conclusions in Sec. 4. 


2 Numerical Method 


2.1 Neural Field Equation without delay 

2.1.1 Time Discretization 

We begin by rewriting equation Q in the form 

c^V(x,t) = I{x,t) -V{x,t) + K{V{x,t)) ( 6 ) 

t e [o,r],x G n C 

where k denotes the nonlinear integral operator defined by 

K{V{x,t))= [ K{\x-y\)S{V{y,t))dy. (7) 

Jn 

We shall first deal with the time discretization in equation (|^, therefore we 
introduce the stepsize ht > 0 and define 

ti = iht, i = 0,...,M, T = htM. 


Moreover, let 

Vi{x) = V{ti,x), Vx G n, i = 0,...,M. 

We shall approximate the partial derivative in time by the backward differ¬ 


ence 


9 T..;. . ^ - 4Vi.i(x) + V-_ 2 (x) 

dt 2ht 


( 8 ) 


which gives a discretization error of the order 0{h^), for sufficiently smooth 
V. By substituting Q into Q we obtain the implicit scheme 

— 4t/j i -|- Ui-2 


2ht 


— li ~ Ui K{Ui), i — 2,..., M, 


(9) 
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where Ui approximates the solution of Q. 

To start this scheme we need to know Uq, which is defined by the initial 
condition Vq, and Ui, which can be obtained by a one-step method, for 
example, the explicit Euler method. 

It is important to analyse the stability of the numerical scheme Q. If 
we denote 

ei = \\Vi-U^\\ 

(the discretization error at time ti), we want to analyse the conditions under 
which, if li = 0, we have Cj < e^-i, if ht is sufficiently small. We can look 
at the scheme as a multistep method for the solution of a system of 
differential equations (the approximate solution Ui , at time ti, is obtained 
from Ui-i and Ui- 2 )- Therefore, we should begin by studying the zero- 
stability. If we multply both sides of Q by 2ht and then ignore the terms 
containing ht, we have an equation of the form 


3Ui — 4:Ui-i + Ui-2 — 0 , 


to which corresponds the characteristic equation 

3//^ — 4/i -|- 1 = 0. 

The roots of this equation are 

2± I 

it is easily seen that \fii^ 2 \ < 1 and = 1 is not a multiple root. Therefore, 
the scheme Q is zero-stable. This means that we have < Cj-i, if ht is 
sufficiently small. 

The above result can be summarized in the form of the following theorem: 

Theorem 2.1. The scheme Q is zero-stable. 

Our next step is to investigate under which conditions equation ([^ has 
a unique solution, so that each step of the iterative process is well defined. 
With this purpose we write this equation in the form 

“ 1 / 3c X (10) 

where 

/,(a;)= + [i, + ^2U-i{x)-^U-2{x)^. xGO (II) 
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Equation @ - 0 is a nonlinear Fredholm integral equation of the second 
kind and we will analyse its solvability using standard results of functional 
analysis. 

In order to apply the Banach hxed point theorem, we define the iterative 
process: 

U^‘'\x) = XK{U^''-^'^) + fi{x) = G{Ur^), xen, n = l,2,... (12) 


were 


A = 


2hf, 


1 ^ 2/if + 3c 


(13) 


If the function G is contractive in a certain closed set X G F, such 
that G{X) C X, then by the Banach fixed point theorem equation (10) has 
a unique solution in X and the sequence , dehned by (12), converges 
to this solution in the norm of F, for any initial guess uf^'^ G X. In our 
case, the solution is by construction the iterate Ui, so it should be close to 
Ui-i and Ui- 2 - Therefore it makes sense to assume that X is a certain set 
containing Ui-i and Ui -2 and to choose = Ui-i. 

To prove that G is contractive in X we need to show that for a certain 
constant L, L < 1, we have 


\\G{V)-G{U)\\f<L\\V-U\\f, ^U,V£X. (14) 


By dehnition 


(||G(B) - G{U)\\Ff = X^ [ \K{x - y)|2|5(y) - SiU)\^dy; (15) 

Jo. 

Using the mean value theorem for integrals, we get 
(||G(U) - G{U)\\Ff < (16) 

U, V t A 


where |0| denotes the area of U. Since we have assumed that S has a 
bounded continuous derivative in M, we can write 

(||G(U) - G{U)\\Ff < |U|A2 (||K||^ 2 (f,.))' Si,, (||U - uyf , (17) 


where 


Smax = max 15'(x) I. 


Finally, we rewrite (17) in the form 

||G(U) - G{U)\\f < X^\\\K\\L2^^2)S^ax\\V - U\\f (18) 
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and we conclude that (14) holds with 




( 19 ) 


Recall that 


A = 


2ht 


< 


2ht 


2ht + 3c '3c 

Then, in order to satisfy L < 1 it is sufficient to require that 


2ht 

3c 


\/W\\\K\\L2(n2)S„ 


< 1 


or equivalently 


ht < 


3c 




( 20 ) 


( 21 ) 


From (21) we conclude that G will be contractive in a certain set X C F 


if we take ht sufficiently small. Moreover, if ht is sufficiently small, we can 
choose a set X such that G{X) C X and {Ui-i, Ui- 2 } C X. Therefore the 


iterative process (12) with = Ut-i will converge to the solution of (10) 


The above construction not only shows that the equation (10) has a 


unique solution in a certain set X, but it also suggests that the iterative 
process (12), starting with = Ut-i can be effectively used to approxi¬ 
mate this solution. Actually, the convergence rate of this process depends 


on the constant L, which as follows from (19) is approximately proportional 


to ht- In other words, the convergence of the process will be faster and faster 
as ht tends to zero. 

The above result can be formulated in the form of the following theorem. 


Theorem 2.2. For each i = 2,3... , if ht satisfies (21) the nonlinear equa¬ 
tion (10) has a unique solution Ui G X, where X C F is a certain closed 


set containing Ut-i and Ui- 2 - Moreover, the iterative process (12) with 
= Ui-iconverges to this solution. 


2.1.2 Space Discretization 


Since the equation (10) in general cannot be solved analytically, we need 


a computational method to compute a numerical approximation of its so¬ 
lution. By other words, we need a space discretization, which will be the 
subject of this subsection. 

For the sake of simplicity, assume that 11 is a rectangle: H = [—1,1] x 
[—1,1]. We now introduce a uniform grid of points {xi,Xj), such that 
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Xi = —1 + ih, i = 0, where h is the discretization step in space. In 
each subinterval we introduce k Gaussian nodes: Xi^s = Xi + ^{1 + 

^s), i = 0,l,...n — 1, where are the roots of the k-ih. degree Legendre 
polynomial, s = We shall denote fl/j the set of all grid points 

(xjs, Xjt)., i,j = 0,n — l,s, t = 1,k. A Gaussian quadrature formula to 
evaluate the integral f{u, v)dudv will have the form 


n—1 n—1 k 


w) = EEEE WsWtf{Xis,Xjt), 


( 22 ) 


i=0 j=0 s=l t=l 


with Ws = ^Ws, where Wg are the standard weights of a Gaussian quadra¬ 
ture formula with k nodes on [—1,1], s = l,...,/c. As it is well-known, a 
quadrature formula of this type has degree 2A: — 1 and therefore, assuming 
that / has at least 2k continuous derivatives on fl, the integration error of 
(22) is of the order of h'^^. Note that the total number of nodes in the space 


discretization is k‘^'n? 


When we introduce the quadrature formula (22) to compute i^{U) we 

Let us denote 


define a finite-dimensional approximation of the operator k 
a vector with entries, where N = nk, such that 


{U 'jisjt ~ U(^Xis, Xjt'), 


then the finite-dimensional approximation of k{U) may be given by 


ni n2 k k 

{K^{U^))mu,lv = ^ ^ X] X] WsWtK{\\{Xmu, Xi^)-{yis, yjt)\\2) S {iU'^)is,jt) ■ 
i=0 1=0 s=l t=l 

(23) 

By replacing n with in equation (10) we obtain the following system of 
nonlinear equations: 


- 


1 


«"([/") = r, 


1 + 1 ^ 


(24) 


where k^{U^) is defined by ( p3| ) and 

(/ )is,jt — f {Xisi Xjt), 

with / defined by (0; in ( |24[ ), for the sake of simplicity, we have omitted 
the index i of 11^. Note that for the the computation of we have to 
evaluate the iterates Ut-i and Ui -2 at all the points of Qh- We denote the 
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vectors resulting from this evaluation by U^_i and U^_ 2 , respectively. We 
conclude that at each time step of our numerical scheme we must solve (24), 
which is a system of nonlinear equations. From this discretization some 
questions arise: 


1. Is the system (24) solvable? 

2. Does the solution of this system converge in some sense to Ui, as 
/i^O? 

3. How can we estimate the error = \\U^ — UA\ ? 


We can investigate the solvability of (24) in the same way as we have studied 
the Fredholm integral equation (10). More precisely, we can introduce the 
iterative procedure 

UhM = + f = (25) 

m = 1, 2,.... As in the case of the Fredholm integral equation, the con¬ 
vergence of the iterative procedure (25) depends on the contractivity of G^. 
Using the same arguments as in subsection 2.1.1, we obtain the following 
inequality, which is a hnite-dimensional analogue of (18): 

n n k k 


\\G\V) - G\U)h < XKrr,a.S'^a.\\y ” Uh E E E E 

2=0 j = 0 t=l 

for U,V eX^ d where 

XfYiax — max |ii"(|| (Xmii, II2) I • 

Since, by the construction of the Gaussian quadrature, 

n n k k 

EEEE'^^'^* = i’ 

2=0 j=0 S=1 t=l 

we get that G^ is Lipschitzian in with the Lipschitz constant 

— -^ATjnax'S'max) 


(26) 


(27) 


(28) 


where A is dehned by (13). Finally we conclude that Gh is contractive if 

3c 


hi < 


2 ATijjax ‘S'max 


(29) 
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Theorem 2.3. For each i = 2, 3... , if S and K are such that Smax o,nd Kmax 
exist and ht satisfies (29), then the nonlinear equation (24) has a unique 
solution G X^, where C is a certain closed set containing Ujfii 
and U^_ 2 - Moreover, the iterative process (25) with = Ujfii converges 

to this solution. 


Since we know that the equation (24) is solvable under certain conditions 
and we even know an iterative scheme to obtain its solution, we should now 
address the problem of the convergence of the solution to C/,, as /i —)• 0. 


With this purpose, we write equation ([1C)|) at each grid point of 

E 


1 


Ufixis^jt) 2/j — ffixis^jt) 1 Xisjt 

1 + T!if 


(30) 


Subtracting from this equation (24), we obtain 

We note that 


i,j = 0, s,t = l,...,k. 

(31) 


(32) 

Substituting (1^ into ([3T| yields 


Ui{xisjt)-Ui,,jt = X{K{Ui{xis,jt))-K. {Ui{is,jt))+K {Ui{is,jt))-K (C/j^jJ), 

(33) 


i,j = 0,..., n — 1, s, t = 1, k. From (33), we obtain 

\\Ui - U’^Woo < X{\\K{Ui) - K^m)\\oc + \\K\U^) - k\U^)\\oo). (34) 

The second term in the last sum can be expressed in terms of \\Ui — U^\\, if 
we take into account that is a Lipschitzian function, that is, there exists 
a certain constant Li such that 


- K^{u^)\\^ < Li\\Ui - u'^\\^,yu^,u^ € X'^, 


where Li is defined by (28). Hence we may rewrite (33) in the form 

\\U^ - u’^Wooii - ALi) < XMUi) - K^{Ui)\U 


(35) 


(36) 


Recall that A is approximately proportional to ht and therefore we have 
ALi < 1 if we choose ht sufficiently small. 
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Finally, in order to evaluate \\K{Ui) —Kp'{Ui))\\oo-, we must remember that 
results from the approximation of the integral k by a Gaussian quadrature 
rule, with stepsize h and k nodes at each subinterval. Therefore, if ?7i, K 
and S are sufficiently smooth, there exists a certain M > 0, which does not 
depend on h, such that 


\K{Ui)-K^{Ui))\\oo<Mh 


2k 


(37) 


Finally, from (36) and (37) we conclude that there exists such a constant M 
that 


\\Ui - U'^Woc < 

The above results lead to the following theorem. 


(38) 


Theorem 2.4. Under the conditions of Theorem 2.3, the unique solution 
of (24) converges to the solution Ui of (10) as /i —)• 0. Moreover, ifUi, 
K and S are sufficiently smooth, the following error estimate holds 


\Ui-U'^\\oc = 0{h^^), as/i^O. 


2.1.3 Computational Implementation 

The above numerical algorithm for the approximate solution of the neural 
field equation in the two-dimensional case was implemented by means of a 
MATLAB code. 

The code has the following structure. After introducing the input data 
(step size in time and in space, initial condition Uq, error tolerance for the 
inner cycle, required number of steps in time) , there is an outer cycle that 
computes each vector , given Ulf_i and Ulf_ 2 , according to the multistep 
method ([^. In order to initialize this cycle, besides Uq, we need C/f, which 
is obtained by the explicit Euler method. More precisely, we compute 


t/f = Uo + -{h-Uo + K^{Uo)). 
c 


(39) 


We recall that at each step in time we must solve the nonlinear system of 


equations (10), which as suggested above is obtained by means of the fixed 


point method, that is, we iterate the scheme (25), until the iterates satisfy 




< e. 


for some given e. This is the inner cycle of our scheme. Typically, in all the 
examples we have computed the number of iterations in the inner cycle is 
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not very high (3-4, in general), confirming that the fixed point method is an 
efficient way of solving the system (10). To start the inner cycle we use an 
initial guess which is obtained from using again the Euler method: 


uh,m = Ui.i + - uti + ( 40 ) 

Note that at each step of the inner cycle it is necessary to compute the 
function at all the grid points. From the computational point of view, 
this means that we must evaluate N'^ times a quadrature rule of the form 
( |22[ ) (with nodes). Of course, this requires a high computational effort 
and the greatest part of the computing time of our algorithm is spent in 
this process. Therefore, we pay special attention to reducing the computa¬ 
tional cost at this stage. In order to improve the efficiency of the numerical 
method, we apply the following technique, proposed in HZj for the solution 
of two-dimensional Fredholm equations. Assuming that the function V is 
sufficiently smooth, we can approximate it by an interpolating polynomial of 
a certain degree. As it is known from the theory of approximation, the best 
approximation of a smooth function by an interpolating polynomial of de¬ 
gree m is obtained if the interpolating points are the roots of the Chebyshev 
polynomial of degree m: 

p^ = cos(——i = (41) 

\ m J 


Our approach for reducing the matrices rank in our method consists in 
replacing the solution V) by its interpolating polynomial at the Chebyshev 
nodes in 11. If V) is sufficiently smooth, this produces a very small error and 
yields a very significant reduction of computational cost. Actually, when 
computing the vector h{Ui) (see formula 0 ) we have only to compute m? 
components, one for each Chebyshev node on [—1,1] x [—1,1]. Choosing m 
much smaller than n, we thus obtain a significant computational advantage. 

The procedure at each iteration is as follows. We compute the matrix 
M such that 


Mij = QiV{p^,p]^,t)), i = 1, = 1, ...,m, 

where Q is the approximation of the integral k, obtained by means of the 
quadrature (22), p^ are the Chebyshev nodes, defined by (41). 

Then we have to perform the matrix multiplication 


A = CMC^, 


(42) 
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where C is the matrix defined by 

Cij = Ci-i{p^), i = j = 

here represents the scaled Chebyshev polynomial of degree k, 

Cfc(x) = cos(A: arcos(x)), A: = 0,1,.. 

with do = Ijy/n, 6k = k = l,...,m — 1. The matrix A contains 

the coefficients of the interpolating polynomial of the solution (expanded in 
terms of scaled Chebyshev polynomials). Finally, in order to obtain the in¬ 
terpolated values of the solution at the Gaussian nodes, we have to compute 

T = P'^AP, (43) 

where P is the transformation matrix, given by 

Pij = Ci-i{xQ)), z = l,...,m, j = l,...,N. 

Here represents each Gaussian node: xq-^ = Xi^s, if j = ik + s. Finally, 
the vector Ui for the next time step (of size ) is obtained by copying T, 
row by row (note that T is a matrix of dimension N x N). 


2.1.4 Complexity Analysis 


As remarked before, it is important to analyse the complexity of the com¬ 
putations, since the computational effort can be signifficantly reduced by 
the application of adequate techniques. In the previous section, we have de¬ 
scribed an algorithm for computing each iterate of the fixed point method, 
which requires applications of the quadrature formula (22). Since this 


quadrature implies evaluations of the integrand function, we have a 
total of m?N‘^ function evaluations. Note that if no polynomial interpola¬ 
tion would be applied, N'^ evaluations of the integrand function would be 
required at each iteration. It is easy to conclude that the number of arith¬ 
metic operations required to apply the quadrature is also proportional to 
rn^N^. 

Then, according to the described algorithm, we must perform the matrix 


multiplication (42). Since the involved matrices have dimension mxm, the 
total number of arithmetic operations is 0{rinA) . Since, by construction, 
m « N, the complexity of this part of the computations is much less than 
the previous one. 
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Finally, we have the matrix product (43) 
trix P has dimensions m x N 


Here the transformation ma- 


as the resulting matrix T has dimensions 
N X N The resulting complexity is therefore 0{mN‘^). 

In conclusion, the number of evaluations of the integrand function in 
each iterate of the hxed point method is N'^m? and the complexity of each 
iteration is 0{N‘^m?). Note that the number of iterations of the fixed point 
method at each time step is typically 2 — 4. 


2.1.5 Error Analysis 

We start by analysing the error resulting from the time discretization. As¬ 
suming that the partial derivatives ^ i = 1)2,3, are continuous on a 

certain domain Q x [0, T], the local discretization error of the approximation, 
given by Q, has the order of 0{h^). 

Concerning the space discretization, the error has two components: one 
resulting from the application of the discretization scheme (24) , and the 


other resulting from the polynomial interpolation. In both cases, the order 
of the approximation depends on the smoothness of the solution. Therefore 
we must choose the degree of the Gaussian quadrature according to the 
smoothness of the mentioned functions. 

The hrst component was analysed in Sec. 2.1.2. According to Theorem 
2.4, if the functions S , K and Ui satisfy certain smoothness conditions, the 


discretization scheme (24) has order 2k in space, where k is the number of 


Gaussian nodes at each subinterval. 

To analyse the interpolation error, we refer to Lemma 3 in mi- Accord¬ 
ing to this Lemma, if the partial derivatives ^ of a certain function 

y3 

f are continuous, with j = 1,2, i = 1, 2,.., s then 


\\f-C^f\\ = 0{m- 


' log^ m) 


(44) 


where Cmf represents the interpolating m-th degree polynomial of / in the 
Ghebyshev nodes. 

In order to obtain an optimal precision of the method, we require that 


the components of the error, given by (38) and (44), are of the same order. 


Hence, to ensure that we can choose m ^ N, the degree of smoothness 
s of the solution should be significantly higher than 2k (otherwise we will 
not obtain a significant reduction of the matrices rank). For example, sup¬ 
pose that m = then for the interpolation error we have 

log^ m = log^(/i“^/^)=—log^ h . Comparing with (38), we 

conclude that for optimal precision we must have s/2 ~ 2k. 
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Hence the described method is specially suitable in the case of a smooth 
solution , so that one can take advantage of the small quadrature error and 
strong rank reduction. 

Summarizing, we have so far shown that the numerical scheme has local 
discretization order 0{h‘f) + 0{h‘^^). Let us now analyse the global error 


For i = 0, we have Eqj = 0,j = 1, For i = 1, we know that Ui is 

computed according to (40). According to the results of section 2.1.2 about 
space discretization, and since the Euler method has local discretization 
error of order 2, we have 


||Ei,,|| = max \Vixj,h) - ([/f),| = 0{ht) + 


(45) 


Concerning the subsequent steps in time, i = 2,3,... from (24) and (10) we 
can conclude that 


4 1 2ht 




As in Section 2.1.2, we can write 

and therefore 


(46) 

(47) 




(48) 


where Li is given by (28). Substituting (48) into (46) we obtain 


which is equivalent to 


\E, 


2h 


1 - ^(1 + Li)) < + 9||E,_2,,|| + 0{h?^), 


i — 2,3, 
(49) 


I 


3c 

provided that 


|((1 + U )< 1 . 


i = 2,3,... 

(50) 

(51) 


In particular, for i = 2, we get 


||i^2,,|| < (l - ^(1 + L)) ' ^||Ei,,|| + ^llEoll + 0(h2^); 


(52) 
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according to (45) we conclude that 


2ht 


\\E2,\\ < (1 - ^(1 + L))-\0{h^^) + 0{ht)) + 0{h^^) = 0{hi) + 0{h^^). 

(53) 

Due to the complexity of our numerical method, it was not possible to obtain 
a closed expression of the global error, for an arbitrary value of i. However, 
according to Theorem 2.1, the multistep scheme is zero-stable, which means 
that it will be stable for a sufficiently small value of /if. The condition (51) 
actually shows us how small /if must be in order to achieve stability. If 
this condition is satisfied, we expect that the scheme will be stable and the 
method will be convergent, with the convergence order 0(/if) + 0{h?^) (the 
same as for the local discretization error). 

We remark that if we used an explicit method, like the Euler method, the 
stability condition on /if would be much more restrictive. In other words, 
the fact that we use an implicit method allows us to use larger step size in 
time, which makes the method more efficient. 

The numerical results presented in Sec. 3 are in agreement with the error 
analysis and confirm the expected convergence orders, for a set of different 
cases. 


2.2 Delay Equation 

We now focus our attention on equation ([^, where the argument of the 
solution inside the integral has a delay r(x, y). This delay takes into account 
the fact that the propagation speed of signals between neurons is finite and 
therefore the post-synaptic potential generated at location x in instant / 
by action potentials arriving from connected neurons at location y actually 
depends on the potential of these neurons at instant t — T{x, y), where t{x, y) 
is the time taken by the signal to come from y to x. Since we assume that 
the propagation speed v is constant and uniform in space, we have 


T{x,y) = 


\y - x\ 


Hence, the delay integro-differential equation that we must solve has the 
form 

= I{x,t)-V{x,t) + f K{\x-y\)S{V{y,t-^-^^—^))dy. (54) 
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Note that in this case the initial conditions satisfied by the solution of our 
problem have the form Q, where 


\y-x 

Tmax = max - 

x,yGQ, V 


The numerical algorithm used to solve equation (54) is essentially the 
same as described in the previous sections. The main difference results from 
the fact that when computing the integral on the right-hand side of (54) at 
instant ti we must use not only the approximate solution at instants and 
ti- 2 , but at all instants k = 1,kmax, where kmax is the integer part 
of Tmax/ht- Note also that the argument ti — may not be a multiple of 
ht- In general let j and 5t be the integer and the fractional part of . In 
this case, we have 

\v — x\ 

ti_j_i <ti- ^ < ti-j 


and 


— { ii 


\y - x\ 




The needed value of the solution V{y,ti — ) is then approximated by 

linear interpolation; 


V 




StUi.j + (1 - 5t)Ui.j.i. 


(55) 


Note that the error analysis that we have carried out in the previous sub¬ 
section may not apply to the delay equation. What we can say in this case, 
assuming that 14 is a smooth function of f, is that the error introduced 
each time we use the approximation formula (55) has the order of 0{h^) 
(the same as the error resulting from the time discretization). However, 
the overall effect of this error in the computations requires a more detailed 
analysis, which is left for a future work. 

Concerning complexity, in the case of the delay equation, each time we 
compute the integrand function, we must compute the delay r(x,y). As 
discussed above, this delay is obtained dividing the distance \\y — x \\2 by 
V. Since this distance is also required to compute the kernel connectivity 
K{\\y — x\\ 2 ), for an effective computation this quantity should be evaluated 
only once and then kept in memory. 
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3 Numerical Results 


3.1 Neural Field Equation without delay 

Here we present the results of some numerical tests we have carried out, 
in order to check the convergence properties of the described method (in 
the case where no delay is considered). Our main purpose is to test exper¬ 
imentally the convergence of the method and measure the error; therefore 
we have chosen some cases where the exact solution is known and do not 
arise from applications. However the form of the connectivity kernels and 
firing rate functions in these examples are close to the ones of neuroscience 
problems. We first check the convergence order in time. With this purpose, 
we consider the following example. 

Example 1. In this example, 

K{\x - y\) = K{xi,X 2 ,yi,y 2 ) = exp {-X{xi - yif - \{x 2 - y 2 f) , 
where A G M+; S{x) = tanh(iTx), a G M+. We set 

/(x, y,t) = — tanh ^aexp ^ y)^ 

where 

6(A,xi,X2) = J^-^J^-^K{xi,X2,yi,y2)dyidy2 = 

= ^ (^Erf(\/A(l - xi)) Erf(\/A(l xi))^ (^Erf(\/A(l - X 2 )) + Erf(\/A(l -h X 2 

where Erf represents the Gaussian error function. 

In this case, it is easy to check that the exact solution is 

V{x,t) = exp(--). 

c 

The initial condition is Vq{x) = 1. 

For the space discretisation, we have used k = 4, that is, 4 Gaussian 
nodes in each subinterval. Since the discretisation error in space must be 
we consider it negligible, compared with the discretisation error in 

time. 

With the following tests, we want to check that the discretisation error 
in time satisfies the condition 

ei = \\Vi-Ui\\=0{hl). 
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t 

ei(O.Ol) 

ei(0.02) 

ei(0.02)/ei(0.01) 

0.02 

6.66E - 5 



0.03 

7.24E - 5 



0.04 

7.46E - 5 

2.66E - 4 

3.57 

0.05 

7.56E - 5 



0.06 

7.61E - 5 

2.91E-4 

3.82 

0.07 

7.65E - 5 



0.08 

7.69E - 5 

3.01E-4 

3.91 

0.09 

7.72E - 5 



0.10 

7.76E - 5 

3.06E - 4 

3.94 


Table 1: Numerical results for Example 1 


The results are displayed in Table 1. We have used two different time 
steps, ht = 0.02 and ht = 0.01, and we have approximated the solution over 
the time interval [0,0.1]. For the space discretisation, we have considered 
N = 24, m = 12. The equation parameters are A = a = c = 1. 

The discretisation errors ei{ht) = jjVi —17j|| are displayed at different mo¬ 
ments ti, for different stepsizes ht. We also present the ratios ei{2ht)/ei{ht), 
which allow us check the convergence order. The ratios are close to 4, which 
conhrms the second order convergence. 

Now, in order to check the convergence of the space discretisation, we 
choose an example, where the time discretisation is exact (does not produce 
any error). 

Example 2. In this example, the functions K and S are the same as in 
example 1. We set 


I{x,y,t) = c + t — tanh((Tf)6(A, x, y). 

As in Example 1, c = 1. In this case, it is easy to check that the exact 
solution is 

V (x, t) = t. 

The initial condition is Vo(x) = 0. The difference operator ([^ is exact 
for linear functions of t, and this is why the scheme in this case does not 
have discretisation error in time. Therefore, the observed errors result from 
the space discretisation. In this case, we are only considering the norm of 
the error at t = 0.1. The derivatives of S and K with respect to the space 
variables depend strongly from the values of A and a, and therefore these 
values should influence the error of the space discretisation. To check this. 
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m 

IV = 12 

iV = 24 

612/624 

= 48 

624/648 

12 

24 

3.11E- 10 

l.llE- 12 

1.03E- 12 

280 

3.997E - 15 
4.413E - 15 

278 

234 


Table 2: Numerical results for Example 2, with A = 1 ,cj = 1. 


m 

II 

to 

4^ 

00 

II 

624/648 

AT = 96 

648/696 

12 

1.62E - 10 

5.52E- 13 

293 

2.36E - 15 

234 

24 

1.69E - 10 

5.33E - 13 

317 

2.22E- 15 

240 


Table 3: Numerical results for Example 2, with A = 5 ,cj = 1. 


we consider 3 different cases: A = 1,(T = 1;A = 1 ,cj = 5; and A = 5, ct = 5, 
which are described in tables 2,3 and 4, respectively. 

Since we are using 4 Gaussian points in each subinterval, we expect 
that the error of the space discretisation is 0{h^) . Therefore, when we 
duplicate the number N of gridpoints, the error should decrease by a factor 
of approximately 2® = 256. 

In order to check the influence of interpolation error, for each N, we 
consider a set of different values of m (interpolation polynomial degree). 

When m increases from 12 to 24, the difference in accuracy is not sig¬ 
nificant. This means that for values of N up to 96 it is enough to consider 
m = 12. When A or ct increase we observe that the errors (for the same 
discretisation step) also increase. This could be expected, since the dis¬ 
cretisation error in space depends on the derivatives ^ and ^Qiy ^ , 

which increase with A and o", respectively. 

Example 3. Einally let us consider an example where the potential 
distribution is not constant, nor in time, neither in space. In this example, 
the function K has the same form as in the previous ones, but the forcing 
function is 

I{xi,X 2 ,t) = - exp /3(A,/U,a;i,X2), 


m 

II 

to 

II 

00 

624/648 

N = 96 

648/696 

12 

7.31E - 10 

2.48E - 12 

295 

9.38E - 15 

264 

24 

7.65E - 10 

2.40E - 12 

319 

8.94E - 15 

268 


Table 4: Numerical results for Example 2, with A = 5, a = 5. 
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ht 

W^ht l|oo 

W^ht oo/ OO 

0.01 

0.005 

0.0025 

7.66E - 5 

1.93E - 5 
4.83E - 6 

3.97 

3.99 


Table 5: Error norms and convergence rates for Example 3, with A = 1, cr = 

1 . 


where 


1 /-i 



0 JO 


exp (-A((xi - yif + {x 2 - y 2 f) - + vl)) dyidy 2 . 


Let us consider a linear firing rate function S{x) = x. If the we set the 
initial condition 

yo{xi,X 2 ) = exp [-y{xl + x|)) , 
we conclude that the exact solution of this problem is 

V(xi,X 2 , t) = exp + ^l)) ■ 


We have computed the numerical solution by our method , over the time 
interval [0,0.05], with stepsize ht = 0.01,0.005 and ht = 0.0025. The pa¬ 
rameters of the space discretisation are m = 12,= 24, therefore the error 
resulting from the space discretisation is in this case negligible, compared 
with the global error. The error analysis for this example is displayed in 
table 5. We see again that the convergence rate is in agreement with the 
theoretical results. 


3.2 Neural Field Equation with Delay 

Example 4. In order to analyse the effect of delay, we now consider equation 
Q with the same data as in Example 3, but with some finite propagation 
speed V. We consider the initial condition Vb(S) = exp (— p(xf -|- x^)), Vx G 

t £ [ Wia2,,0]. 

In Fig. 1 some graphs of the solution are displayed. On the left-hand 
side, one can see the plots of the solution at t = 0.5, t = 1, t = 1.5, and 
t = 2.0, for the non-delay case. On the right-hand side the corresponding 
plots are depicted, but for the delayed equation, when the propagation speed 
is u = 1. The values of the remaining parameters are c=l,/U = l,A = l. 
In both cases the stepsize in time is ht = 0.1, and the parameters of the 
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space discretisation are m = 12, = 24. From this figure it is clear that 

as an effect of the delay, the decay of the solution in the case of the delayed 
equation is much slower. 


4 Conclusions 

We have described and analysed a new numerical algorithm for computing 
approximate solutions of the two-dimensional neural field equations with de¬ 
lay. The stability, convergence and complexity of this algorithm have been 
analysed and a set of numerical examples has been presented. The numerical 
experiments are in agreement with the theoretical results and confirm that 
this algorithm can be successfully used for the solution of problems in Neu¬ 
roscience and Robotics. The main advantages of the described algorithm are 
its stability and accuracy, which is illustrated by the presented examples. 
Moreover, due to the use of a rank reduction technique, it is efficient when 
dealing with the two-dimensional case. 
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